source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

fdx1.path <- list.dirs('mission/fdx1', recursive = F)

fdx1.mex <- fdx1.path |>
  str_extract('KO|WT') |>
  set_names(fdx1.path, nm = _) |>
  Read10X()

fdx1.mex[1:5,1:5]

sobj <- fdx1.mex |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

sobj %<>% PercentageFeatureSet('^mt-', col.name = 'mito.ratio')

sobj |> VlnPlot('mito.ratio', pt.size = 0) +
  geom_hline(yintercept = 20, linetype = 'dashed')

sobj |> ggplot(aes(orig.ident, mito.ratio)) +
  geom_violin(draw_quantiles = 0.95)

sobj.clean <- sobj |> filter(mito.ratio < 20)

sobj.clean %<>% quick_process_seurat()

mmur <- celldex::MouseRNAseqData()

sobj.clean %<>% mark_cell_type_singler(mmur, new_label = 'mmur.main')

sobj.clean |> DimPlot(group.by = 'mmur.main', split.by = 'orig.ident')

gut.marker <- list(
  stem.cells = c('Lgr5','Hmgcs2','Cdx2','Axin2','Ascl2','Fabp2','Ugt2b5'),
  DCS = c('Reg4','Tff3','Kit','Mptx1','Sval1','Spink1','Best2','Fhl1','Clec2h'),
  goblet.cells = c('Spdef','Ccl9','B3gnt6','Casc4','Smim14','Plet1','Cgref1'),
  TA.cells = c('Mki67','Top2a','Cenpf','Pclaf','Cenpa','Cdk1','Ube2c'),
  Colonocytes = c('Alpi','Slc6a14','Wfdc2','Sult1b1','Maoa','Fgfbp1','Mep1a'),
  Enteroendocrine = c('Chga','Chgb','Neurod1','Neurog3','Gfra3'),
  Tuft.cells = c('Col1a1','Col1a2','Col3a1','Igfbp7','Serpinh1','Dclk1'))

sobj.clean |> 
  filter(orig.ident == 'WT') |>
  DotPlot(gut.marker,cluster.idents = T) +
  RotatedAxis()

sobj.clean %<>%
  mutate(cell.subset = case_match(as.integer(seurat_clusters),
                                  7 ~ 'DCS.cells',
                                  6 ~ 'Goblet.cells',
                                  16 ~ 'Enteroendocrine',
                                  c(3,11) ~ 'TA.cells',
                                  c(1,9,10,4) ~ 'Stem.cells',
                                  c(8,5) ~ 'Colonocytes',
                                  .default = mmur.main))

sobj.clean |> DimPlot(group.by = 'cell.subset', cols = DiscretePalette(36),
                      label = T)

clus3 <- FindMarkers(sobj.clean, ident.1 = 3, logfc.threshold = 1) |>
  as_tibble(rownames = 'gene')

clus3 |> filter(avg_log2FC > 0)

# rds checkpoint ------
sobj.clean |> write_rds('mission/fdx1/gutepi.rds')

sobj.clean <- read_rds('mission/fdx1/gutepi.rds')

sobj.clean |>
  DotPlot(features = gut.marker, group.by = 'cell.subset')

Idents(sobj.clean) <- 'cell.subset'

my.marker <- sobj.clean |>
  FindAllMarkers(logfc.threshold = 1, only.pos = T)

marker5 <- my.marker |>
  as_tibble() |>
  filter(p_val_adj < .01, str_detect(gene, 'Gm|ENSMU|Rik', negate = T)) |>
  distinct(gene, .keep_all = T) |>
  group_by(cluster) |>
  slice_min(p_val_adj, n = 5, with_ties = F) |>
  arrange(as.character(cluster))

sobj.clean |>
  ScaleData(features = marker5$gene) |>
  DoHeatmap(features = marker5$gene, group.by = 'cell.subset',
            label = F, group.colors = DiscretePalette(12)) +
  scale_fill_viridis_c()

# cell type fraction -------
type.change <- sobj.clean |>
  as_tibble() |>
  discov_frac_change(orig.ident, cell.subset, KO, WT)

type.change |>
  ggplot(aes(subtype, log2fc_frac, fill = type)) +
  geom_col() +
  theme_pubr(x.text.angle = 45,base_size = 14) +
  scale_fill_manual(values = c('blue','red','grey')) +
  labs(title = 'Fraction changes of cell types in KO',
       y = 'Log2fc of fraction')

type.order <- type.change |>
  mutate(cell.subset = as.character(subtype),
         log2fc_frac, .keep = 'none')

sobj.clean |>
  as_tibble() |>
  left_join(type.order) |>
  mutate(group = fct_relevel(orig.ident, 'KO'),
         cell.subset = fct_reorder(cell.subset, log2fc_frac)) |>
  ggplot(aes(y = group, fill = cell.subset)) +
  geom_bar(position = 'fill') +
  scale_fill_brewer(palette = 'Paired') +
  theme_pubr(legend = 'right') +
  labs(x = 'Fraction')

type.conf <- sobj.clean |>
  calc_frac_conf_on_grouped_count(orig.ident, cell.subset)

type.conf |>
  left_join(type.order) |>
  mutate(cell.subset = fct_reorder(cell.subset, log2fc_frac),
         group = fct_relevel(orig.ident, 'WT')) |>
  ggplot(aes(group, fraction, ymax = conf.high, ymin = conf.low,
             fill = group)) +
  geom_col() +
  geom_errorbar(width = .5) +
  facet_wrap(~cell.subset, scales = 'free_y') +
  scale_fill_hue(direction = -1) +
  theme_pubr(legend = 'right') +
  labs(title = 'Fraction changes of cell types in intestine epithelium')

# DEGs ------
## Fdx1 -----
sobj.clean |>
  mutate(subgroup = str_c(cell.subset, orig.ident, sep = '_')) |>
  DotPlot('Fdx1', group.by = 'subgroup')

sobj.clean |>
  DotPlot('Fdx1', split.by = 'orig.ident')

sobj.clean |>
  mutate(orig.ident = fct_relevel(orig.ident, 'WT')) |>
  VlnPlot('Fdx1', group.by = 'cell.subset', split.by = 'orig.ident',
          cols = c('blue','red'), pt.size = 0)

sobj.clean |>
  get_abundance_sc_long('Fdx1') |>
  left_join(x = sobj.clean, y = _) |>
  mutate(group = fct_relevel(orig.ident, 'WT')) |>
  ggplot(aes(cell.subset, .abundance_RNA, fill = group, color = group)) +
  geom_violin(scale = 'width') +
  theme_pubr(x.text.angle = 45) +
  scale_fill_hue(direction = -1) +
  scale_color_hue(direction = -1) +
  stat_summary(fun = 'ExpMean', geom = 'crossbar') +
  labs(title = 'Fdx1 RNA expression', y = 'Normalized expression')

## cell.subset -----
subset.list <- sobj.clean$cell.subset |> unique()

subset.list %<>%
  set_names(nm = subset.list)

subset.kovwt.deg <- subset.list |>
  map(\(x)FindMarkers(sobj.clean, subset.ident = x, min.cells.group = 1,
                    group.by = 'orig.ident', ident.1 = 'KO',) |>
        as_tibble(rownames = 'gene'), .progress = T) |>
  list_rbind(names_to = 'cell.subset')

subset.kovwt.deg |>
  filter(p_val_adj < .05) |>
  ggplot(aes(cell.subset, fill = cell.subset)) +
  geom_bar() +
  scale_fill_brewer(palette = 'Paired') +
  theme_pubr(legend = 'right', x.text.angle = 45) +
  ggtitle('Number of significant differential expressed genes in KO vs WT')

volc.list <- subset.list |>
  map(\(x)subset.kovwt.deg |>
        filter(cell.subset == x) |>
        plot_bill_volc(group1 = 'KO', group2 = 'WT') +
        ggtitle(x))

volc.list[[4]]

subset.kovwt.deg |> write_csv('mission/fdx1/subset.kovwt.deg.csv')

## pathway -----
library(clusterProfiler)

### ORA -------
up.ora.go.list <- subset.list |>
  map(\(x)subset.kovwt.deg |>
        filter(cell.subset == x, p_val_adj < .05, avg_log2FC > 0) |>
        pull(gene) |>
        enrichGO(OrgDb = 'org.Mm.eg.db',
                 keyType = 'SYMBOL',
                 ont = 'ALL',
                 minGSSize = 3, readable = T) |>
        pluck('result'), .progress = T)

up.ora.go.list <- up.ora.go.list |>
  map(as_tibble, rownames = 'id') |>
  list_rbind(names_to = 'cell.subset')

up.ora.go.list |>
  filter(p.adjust < .05, qvalue < .05, Count > 2) |>
  write_csv('mission/fdx1/ko.ora.go.list.csv')

up.ora.go.list |>
  filter(str_detect(cell.subset, 'Colono')) |>
  slice_max(Count, n = 5, by = ONTOLOGY) |>
  ggplot(aes(y = Description, x = Count, fill = qvalue)) +
  geom_col() +
  facet_wrap(~ONTOLOGY, scales = 'free_y', ncol = 1,strip.position = 'left') +
  scale_fill_gradient(low = 'red', high = 'black') +
  theme_pubr(legend = 'right') +
  labs(title = 'Upregulated GO pathway in Colonocytes in KO')

highlighted.stem <- c('Cd177','Duox2','Plet1','Duoxa2','Mcpt1','Mcpt2','Igkc')
colono.high <- c('Wfdc18','Clca4b','Cd177','Ly6e','Clca3b','Khdc1a','Ly6d','Rbp1','Prss22')
colono.low <- c('Klk1','Sult1a1','Fosb','Cyp2c68','Selenbp1','Mptx1','Zg16','Sycn','Cyp2c69')

up.ora.go.list |>
  filter(Count > 2) |>
  separate_longer_delim(geneID, delim = '/') |>
  filter(geneID %in% colono.high,
         str_detect(cell.subset, 'Colon')) |>
  dplyr::count(Description,sort = T)

up.ora.go.list |>
  filter(str_detect(cell.subset, 'DCS'), Count > 1) |>
  plot_enrichment(n = 8) +
  ggtitle('Upregulated GO pathway in DCS cells in KO')
  
sobj.clean |>
  filter(str_detect(cell.subset, 'TA')) |>
  VlnPlot(features = c('Gpx2','Duoxa2','Reg3b'),
            group.by = 'orig.ident')

sobj.clean |>
  filter(str_detect(cell.subset, 'TA')) |>
  get_abundance_sc_long(features = c('Gpx2','Duoxa2','Reg3b')) |>
  mutate(group = str_extract(cell, 'KO|WT') |> fct_relevel('WT')) |>
  ggplot(aes(group, .abundance_RNA, fill = group)) +
  geom_violin() +
  stat_summary(fun = 'ExpMean', geom = 'crossbar', width = .5, fill = 'black') +
  facet_wrap(~.feature) +
  theme_pubr(legend = 'right') +
  scale_fill_hue(direction = -1) +
  labs(title = 'Upregulated "Regulation of inflammatory response"\nGO pathway in KO TA cells')

  
down.ora.go.list <- subset.list |>
  map(\(x)subset.kovwt.deg |>
        filter(cell.subset == x, p_val_adj < .05, avg_log2FC < 0) |>
        pull(gene) |>
        enrichGO(OrgDb = 'org.Mm.eg.db',
                 keyType = 'SYMBOL',
                 ont = 'ALL',
                 minGSSize = 3, readable = T) |>
        pluck('result'), .progress = T)

down.ora.go.list <- down.ora.go.list |>
  map(as_tibble, rownames = 'id') |>
  list_rbind(names_to = 'cell.subset') |>
  filter(p.adjust < .05, qvalue < .05, Count > 2) |>
  write_csv('mission/fdx1/wt.ora.go.list.csv')
  
down.ora.go.list |>
  filter(str_detect(cell.subset, 'Stem')) |>
  slice_max(Count, n = 5, by = ONTOLOGY, with_ties = F) |>
  ggplot(aes(y = str_wrap(Description, width = 40), x = Count, fill = qvalue)) +
  geom_col() +
  facet_wrap(~ONTOLOGY, scales = 'free_y', ncol = 1,strip.position = 'left') +
  scale_fill_gradient(low = 'blue', high = 'black') +
  theme_pubr(legend = 'right') +
  labs(title = 'Downregulated GO pathway in Stem cells in KO',
       y = 'Description')

down.stem <- c('S100g','Klk1','Mptx1','Oat','Slc30a10','Mtus2','Capsl','Cyp4b1','Rims2')

down.ora.go.list |>
  filter(Count > 2) |>
  separate_longer_delim(geneID, delim = '/') |>
  filter(geneID %in% down.stem,
         str_detect(cell.subset, 'Stem')) |>
  dplyr::count(Description,sort = T)

up.ora.kegg.list <- subset.list |>
  map(\(x)subset.kovwt.deg |>
        filter(cell.subset == x, p_val_adj < .05, avg_log2FC > 0) |>
        pull(gene) |>
        bitr(fromType = 'SYMBOL', toType = 'ENTREZID',
             OrgDb = 'org.Mm.eg.db') |>
        pull(ENTREZID) |>
        enrichKEGG(organism = 'mmu', keyType = 'kegg') |>
        pluck('result') |>
        as_tibble(rownames = 'id'), .progress = T)

up.ora.kegg.list <- up.ora.kegg.list |>
  map(as_tibble, rownames = 'id') |>
  list_rbind(names_to = 'cell.subset') |>
  write_csv('mission/fdx1/ko.ora.kegg.list.csv')

up.ora.kegg.list |>
  mutate(Description = str_remove(Description, ' - .+')) |>
  filter(str_detect(cell.subset, 'Colon')) |>
  plot_enrichment(n = 15)

stem.down.ora.go <- subset.kovwt.deg |>
  filter(cell.subset == 'Stem.cells', p_val_adj < .05, avg_log2FC > 0) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'ALL',
           minGSSize = 3, readable = T)

stem.up.ora.go |> pluck('result') |> as_tibble()

# heatmap of SDEG in stem/colono --------
subset.kovwt.deg <- read_csv('mission/fdx1/subset.kovwt.deg.csv')

sdeg.stemcolo <- subset.kovwt.deg |>
  filter(str_detect(cell.subset, 'Stem|Colon'), p_val_adj < .05,
         abs(avg_log2FC) > 1) |>
  summarise(n = n(), .by = 'gene') |>
  filter(n > 1) |>
  pull(gene)

# IFN pathways ------
ifn.list <- sobj.clean |> rownames() |> str_subset('Ifna|Infb')

sobj.clean |>
  mutate(subgroup = str_c(cell.subset, '_', orig.ident)) |>
  DotPlot(c('Ifna2','Zbp1'), group.by = 'subgroup')

# mast cells ------
mast.marker <- c('Kit','Cma1','Cpa3','CTSG','HDC','MCPT1','MCPT4','MCPT2','TPSB2') |>
  str_to_title()

sobj.clean |> DotPlot(mast.marker, group.by = 'cell.subset') +
  labs(title = 'Mast cell marker gene', x = 'Gene', y = 'Cell type')

sobj.clean |> DimPlot(cols = DiscretePalette(36))

sobj.clean |> DotPlot(mast.marker, group.by = 'seurat_clusters')

sobj.clean |> FeaturePlot(c('Kit','Cpa3','Mcpt1','Mcpt2'))

sobj.clean |>
  filter(str_detect(cell.subset, 'Colonocyte'))

sobj.clean |>
  get_abundance_sc_long('Mcpt1') |>
  left_join(sobj.clean, y = _) |>
  ggplot(aes(cell.subset, .abundance_RNA)) +
  stat_summary(geom = 'col') +
  theme_pubr(x.text.angle = 45) +
  labs(title = 'Average expression of Mcpt1') +
  facet_wrap(~orig.ident)

# cell death pathway in epi --------
autophagy <- c('Ulk2','Atg14','Tex264','Atg101','Atg7',
               'Becn1','Wipi1','Map1lc3b')

pyroptosis <- c('Aim2','Casp1','Casp4','Nlrp3','Nlrp1b','Gsdmd','Mefv')

ferroptosis <- c('Gss','Slc40a1','Gpx4','Ascl4','Pcbp1','Pcbp2')

necroptosis <- c('Mlkl','Cyld','Ripk3','Ripk1','Zbp1')

sobj.clean$cell.subset |> table()

epi.osis <- sobj.clean |>
  filter(str_detect(cell.subset, 'Colo|DCS|Entero|Gob|Ste|TA')) |>
  get_abundance_sc_long(c(autophagy, pyroptosis, ferroptosis, necroptosis))

## average expr / scaled average ---------
epi.osis.scal <- epi.osis |>
  mutate(group = str_extract(cell, 'KO|WT'),
         foo = ifelse(str_ends(cell, 'A|G'), '1','2'),
         subgroup = str_c(group, foo)) |>
  summarise(avg.expr = log1p(ExpMean(.abundance_RNA)),
            .by = c(.feature, group)) |>
  mutate(scaled.expr = scale(avg.expr)[,1], .by = .feature)

epi.osis.scal |>
  write_csv('mission/fdx1/cell.death.average.csv')

epi.osis |>
  distinct(cell) |>
  dplyr::count(str_detect(cell, 'KO'), str_detect(cell, 'A$|G$'))

epi.osis.scal |>
  mutate(pathway = case_match(.feature,
                              autophagy ~ 'autophagy',
                              ferroptosis ~ 'ferroptosis',
                              pyroptosis ~ 'pyroptosis',
                              .default = 'necroptosis'),
         subgroup = fct_relevel(subgroup, 'WT1','WT2')) |>
  ggplot(aes(subgroup, .feature, fill = scaled.expr)) +
  geom_tile(color = 'black') +
  scale_fill_distiller(palette = 'RdYlBu') +
  facet_wrap(~pathway, scales = 'free', ncol = 4) +
  theme_pubr(legend = 'right') +
  labs(title = 'Cell death pathway in Fdx1-cKO intestine epithilium',
       y = 'gene')

scale(1:5)[,1]

deaths <- c('autophagy', 'ferroptosis', 'pyroptosis', 'necroptosis')

death.plots <- deaths |>
  map(\(x)epi.osis.scal |>
        mutate(subgroup = fct_relevel(subgroup, 'WT1','WT2')) |>
        filter(.feature %in% get(x)) |>
        ggplot(aes(subgroup, .feature, fill = scaled.expr)) +
        geom_tile(color = 'black') +
        scale_fill_distiller(palette = 'RdYlBu') +
        theme_pubr(legend = 'right') +
        labs(title = x, y = 'gene'))

patchwork::wrap_plots(death.plots,ncol = 4)

death.plots.avg <- deaths |>
  map(\(x)epi.osis.scal |>
        mutate(group = fct_relevel(group, 'WT')) |>
        filter(.feature %in% get(x)) |>
        ggplot(aes(group, .feature, fill = avg.expr)) +
        geom_tile(color = 'black') +
        scale_fill_distiller(palette = 'RdYlBu') +
        theme_pubr(legend = 'right') +
        labs(title = x, y = 'gene'))

death.plots.avg[[1]]

patchwork::wrap_plots(death.plots.avg,ncol = 4)

## doheatmap ------
scheat.death <- deaths |>
  map(\(x)sobj.clean |>
        ScaleData(features = get(x)) |>
        DoHeatmap(features = get(x),group.by = 'orig.ident') +
        scale_fill_distiller(palette = 'RdYlBu') +
        ggtitle(x))

patchwork::wrap_plots(scheat.death)

## fold change ------
epi.logfc <- sobj.clean |>
  filter(str_detect(cell.subset, 'Colo|DCS|Entero|Gob|Ste|TA')) |>
  FindMarkers(group.by = 'orig.ident', ident.1 = 'KO', logfc.threshold = 0,
              min.pct = 0,
              features = c(autophagy, pyroptosis, ferroptosis, necroptosis))

### heatmap
tidy.logfc <- epi.logfc |>
  as_tibble(rownames = 'gene') |>
  filter(gene %in% c(autophagy, pyroptosis, ferroptosis, necroptosis)) |>
  mutate(pathway = case_match(gene,
                              autophagy ~ 'autophagy',
                              ferroptosis ~ 'ferroptosis',
                              pyroptosis ~ 'pyroptosis',
                              .default = 'necroptosis'),
         group = 'KO')

tidy.logfc |>
  mutate(group = 'WT', avg_log2FC = 0) |>
  bind_rows(tidy.logfc) |>
  mutate(group = fct_relevel(group, 'WT')) |>
  ggplot(aes(group, gene, fill = avg_log2FC)) +
  geom_tile(color = 'black') +
  facet_wrap(~pathway, scales = 'free_y', ncol = 4) +
  scale_fill_distiller(palette = 'RdYlBu', values = c(0,.25,.5,.6,1)) +
  theme_pubr(legend = 'right')

tidy.logfc |>
  mutate(group = 'WT', avg_log2FC = 0) |>
  bind_rows(tidy.logfc) |>
  pivot_wider(names_from = group, values_from = avg_log2FC) |>
  select(-c(2:4)) |>
  write_csv('mission/fdx1/cell.death.epithelium.csv')

### barplot
epi.logfc |>
  as_tibble(rownames = 'gene') |>
  filter(gene %in% c(autophagy, pyroptosis, ferroptosis, necroptosis)) |>
  mutate(pathway = case_match(gene,
                              autophagy ~ 'autophagy',
                              ferroptosis ~ 'ferroptosis',
                              pyroptosis ~ 'pyroptosis',
                              .default = 'necroptosis'),
         type = case_when(avg_log2FC > 0 & p_val_adj < .05 ~ 'Up',
                          avg_log2FC < 0 & p_val_adj < .05 ~ 'Down',
                          .default = 'NS')) |>
  ggplot(aes(gene, avg_log2FC, fill = type)) +
  geom_col() +
  facet_wrap(~pathway, scales = 'free')

epi.logfc |>
  as_tibble(rownames = 'gene') |>
  filter(gene %in% c(necroptosis)) 
